options(width=108)
out=F
# The list of subjects, the order of conditions, and the thresholds are derived from Subjects.xlsx
library(xlsx)
library(ggplot2)
library(plyr)
library(reshape)
library(matrixStats)
library(gridExtra)
library(lme4)
library(lmerTest)
library(BayesFactor)
#library(splines)
db <- '/home/egor/Dropbox/' # on Linux
db <- '/Users/Egor/Dropbox/' # Windows
# db <- '~/Dropbox/' # on Mac
source(paste(db, 'Prog/R/myFunctions/pvalfn.R', sep=''))
# theme for plotting:
alpha <- .7
w <- .56 # proportion width in group plots
themefy <- function(p) {
p <- p + theme_bw() +
theme(panel.grid.minor.x=element_blank(), panel.grid.minor.y=element_blank(),
axis.text=element_text(size=8), axis.title=element_text(size=9),
legend.text=element_text(size=8), legend.title=element_text(size=9),
legend.key = element_blank(), legend.margin=margin(t=-.04, unit='in'),
legend.background = element_rect(fill='transparent'),
plot.title=element_text(face='bold'))
}
cc <- c('#e31a1c','#fdbf6f','#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99')
xLab <- 'Target Peak Time (s)'
yLab <- 'Log Contrast Threshold'
# colLab <- expression(paste('\nTarget\nVelocity (', degree, '/s)', sep=''))
colLab <- expression(paste('Target Eccentricity (', degree, ')', sep=''))
# colLabType <- 'Mask Type'
yLim <- c(-1.5,-0.5)
dodge <- position_dodge(width=0.1)
sumFn <- function(ss, subjStr='subj', xStr='targTpeak', grpStr='targEcc'){
sumSubj <- ddply(ss, c(subjStr, xStr, grpStr), summarise,
mnS=mean(threshMean), se=sd(threshMean)/sqrt(length(threshMean))) #, .drop=F)
sumSubjMn <- ddply(ss, c(subjStr), summarise, mnStot=mean(threshMean)) # total mean across conditions per subj
sumSubj <- merge(sumSubj, sumSubjMn)
sumSubj$normS <- - sumSubj$mnS / sumSubj$mnStot # normalized subject mean
sumSubj$seNorm <- NA
# sumSubj$mnS[is.na(sumSubj$mnS)] <- 0
sumGrp <- ddply(sumSubj, c(xStr, grpStr), summarise,
mn=mean(mnS), se=sd(mnS)/sqrt(length(mnS)),
norm=mean(normS), seNorm=sd(normS)/sqrt(length(normS)))
sumGrp$subj <- 'average'
sumSubj <- rename(sumSubj, c(mnS='mn',normS='norm'))
sumComb <- rbind(sumGrp, subset(sumSubj, select=-mnStot))
sumComb$se[is.na(sumComb$se)] <- 0
sumComb
}
plotAve <- function(pss, subjStr='subj', xStr='targTpeak', grpStr='targEcc',
xlab=xLab, ylab=yLab, collab=colLab, yStr='mn', seStr='se'){
pss$yMin <- pss[,yStr] - pss[,seStr]
pss$yMax <- pss[,yStr] + pss[,seStr]
pss[,grpStr] <- factor(pss[,grpStr])
p <- ggplot(pss, aes_string(x=xStr, y=yStr, colour=grpStr, group=grpStr,
ymin='yMin', ymax='yMax')) +
geom_point(position=dodge, size=1, alpha=alpha) + geom_line(position=dodge, alpha=alpha) +
scale_x_continuous(breaks=c(.5,1,1.5), labels=c('0.5','1','1.5')) +
geom_linerange(position=dodge, show.legend=F, alpha=alpha) +
labs(x=xlab, y=ylab, colour=collab) + #ylim(0,1) +
guides(colour=guide_legend(keyheight=.3, default.unit='inch'))
p <- themefy(p)
}
plotIndiv <- function(pss, subjStr='subj', xStr='targTpeak', grpStr='targEcc',
xlab=xLab, ylab=yLab, collab=colLab, yStr='mn', seStr='se'){
pss$yMin <- pss[,yStr] - pss[,seStr]
pss$yMax <- pss[,yStr] + pss[,seStr]
pss[,grpStr] <- factor(pss[,grpStr])
p <- ggplot(pss, aes_string(x=xStr, y=yStr, colour=grpStr, group=grpStr,
ymin='yMin', ymax='yMax')) +
facet_wrap( ~ subj, ncol=4) +
geom_point(position=dodge, size=1, alpha=alpha) + geom_line(position=dodge, alpha=alpha) +
geom_linerange(position=dodge, show.legend=F, alpha=alpha) +
scale_x_continuous(breaks=c(.5,1,1.5), labels=c('0.5','1','1.5')) +
# scale_y_continuous(breaks=c(0,.25,.5,.75,1), labels=c('0','','0.5','','1'), limits=c(0,1)) +
labs(x=xlab, y=ylab, colour=collab) +
guides(colour=guide_legend(keyheight=.3, default.unit='inch'))
p <- themefy(p)
}
allDataDir <- paste(db,'Projects/mc/data_cfs/cfs',sep='')
dataDirs <- dir(allDataDir)
dataDirs <- dataDirs[grep('cfs',dataDirs)]
colsOfInt <- c('participant', 'dom', 'session', 'mcBv', 'targTpeak', 'targXoff2',
'targV', 'stairStart', 'meanRev6')
df <- data.frame()
dfRevs <- data.frame()
dfIntn <- data.frame()
for(curDir in dataDirs){
print(curDir)
curDf <- read.csv(paste(allDataDir,'/',curDir,'/',curDir,'.csv',sep=''))
curDf <- curDf[,colsOfInt]
curDf$cond <- substr(curDir,19,22)
subjStairs <- dir(paste(allDataDir,'/',curDir,'/',sep=''))
subjStairs <- subjStairs[grep('.tsv',subjStairs)]
for(curStairFN in subjStairs){
curStair <- paste(allDataDir,'/',curDir,'/',curStairFN,sep='')
curRevs <- read.table(curStair, skip=1, nrows=1)
curIntn <- read.table(curStair, skip=4, nrows=2)
curInfo <- readLines(curStair)
curInfoCols <- data.frame(subj=curDf$participant[1], dom=curDf$dom[1],
sess=curDf$session[1], stairStart=curInfo[33],
mcBv=curInfo[41], targTpeak=curInfo[43],
targEcc=curInfo[37], targV=curInfo[31])
curDfRevs <- cbind(curInfoCols, curRevs[,2:11])
nTrials <- ncol(curIntn)-1
curDfIntn <- curInfoCols[rep(seq_len(nrow(curInfoCols)), each=nTrials),]
curDfIntn$trial <- seq(1,nTrials)
curDfIntn$intn <- as.numeric(curIntn[1,2:(nTrials+1)])
curDfIntn$resp <- as.numeric(curIntn[2,2:(nTrials+1)])
rownames(curDfIntn) <- NULL
dfRevs <- rbind(dfRevs, curDfRevs)
dfIntn <- rbind(dfIntn, curDfIntn)
}
df <- rbind(df, curDf)
}
## [1] "mc2_cfs_cent_p12_s2_2017-11-02_1537"
## [1] "mc2_cfs_cent_p13_s2_2017-11-01_1030"
## [1] "mc2_cfs_cent_p17_s1_2017-10-26_1039"
## [1] "mc2_cfs_cent_p18_s1_2017-10-26_1443"
## [1] "mc2_cfs_dyna_p10_s2_2017-10-27_1511"
## [1] "mc2_cfs_dyna_p11_s1_2017-11-02_1439"
## [1] "mc2_cfs_dyna_p14_s1_2017-10-27_1235"
## [1] "mc2_cfs_dyna_p4_s2_2017-10-26_1342"
## [1] "mc2_cfs_peri_p12_s1_2017-10-26_1539"
## [1] "mc2_cfs_peri_p13_s1_2017-10-31_1032"
## [1] "mc2_cfs_peri_p17_s2_2017-11-02_1018"
## [1] "mc2_cfs_peri_p18_s2_2017-11-02_1327"
## [1] "mc2_cfs_stat_p10_s1_2017-10-27_1030"
## [1] "mc2_cfs_stat_p11_s2_2017-11-02_1630"
## [1] "mc2_cfs_stat_p14_s2_2017-10-27_1319"
## [1] "mc2_cfs_stat_p4_s1_2017-10-26_1234"
ds <- rename(df, c(participant='subj', session='sess', meanRev6='thresh',
mcBv='maskV', targXoff2='targEcc'))
ds$targEcc <- round(ds$targEcc / 35,1)
ds$targV <- round(ds$targV / 3.5,1)
ds$maskType <- ''
ds$maskType[ds$maskV==0] <- 'static'
ds$maskType[ds$maskV==10] <- 'standard'
ds$maskType[ds$maskV==60] <- 'fast'
# ds$maskV <- round(ds$maskV * 60 / 35, 1)
# ds$maskV[ds$maskV<0.05] <- 0
ds$condSplit <- ''
ds$condSplit[ds$cond=='stat' | ds$cond=='dyna'] <- 'v'
ds$condSplit[ds$cond=='cent' | ds$cond=='peri'] <- 'ecc'
head(ds)
## subj dom sess maskV targTpeak targEcc targV stairStart thresh cond maskType condSplit
## 1 12 1 2 60 1.5 0.8 4.6 0 -1.378333 2_20 fast
## 2 12 1 2 0 1.5 0.8 4.6 -2 -1.375000 2_20 static
## 3 12 1 2 60 0.5 0.8 0.0 -2 -1.263333 2_20 fast
## 4 12 1 2 60 1.0 0.8 4.6 -2 -1.341667 2_20 fast
## 5 12 1 2 10 1.5 0.8 0.0 0 -1.363333 2_20 standard
## 6 12 1 2 0 1.5 0.8 0.0 0 -1.636667 2_20 static
thresh <- ddply(ds, .(subj,dom,sess,maskV,maskType,targTpeak,targEcc,targV,
cond,condSplit), summarise, threshMean = mean(thresh),
threshSD = sd(thresh))
head(thresh)
## subj dom sess maskV maskType targTpeak targEcc targV cond condSplit threshMean threshSD
## 1 4 1 1 0 static 0.5 0.8 0 _201 -1.565833 0.11667262
## 2 4 1 1 0 static 0.5 2.9 0 _201 -1.122500 0.14024284
## 3 4 1 1 0 static 1.0 0.8 0 _201 -1.219167 0.02474874
## 4 4 1 1 0 static 1.0 2.9 0 _201 -0.910000 0.30876996
## 5 4 1 1 0 static 1.5 0.8 0 _201 -1.342500 0.08367430
## 6 4 1 1 0 static 1.5 2.9 0 _201 -1.042500 0.13552880
ss <- sumFn(thresh[(thresh$maskType=='static' & thresh$targV==0),])
ssAve <- ss[ss$subj=='average',]
ssIndiv <- ss[ss$subj!='average',]
pTargStatMaskStat <- themefy(plotAve(ssAve)) + ylim(yLim)
(plotIndiv(ssIndiv) + ylim(-2,0))
iss <- dfIntn[dfIntn$mcBv==0 & dfIntn$targV==0,]
(p <- ggplot(iss, aes(x=trial, y=intn, colour=targEcc, linetype=stairStart)) +
facet_wrap(subj ~ targTpeak, ncol=6) + geom_point() + geom_line() + ylim(-2,0))
ss <- sumFn(thresh[thresh$maskType=='static' & thresh$targV>0,])
ssAve <- ss[ss$subj=='average',]
ssIndiv <- ss[ss$subj!='average',]
pTargDynaMaskStat <- themefy(plotAve(ssAve)) + ylim(yLim)
(plotIndiv(ssIndiv) + ylim(-2,0))
iss <- dfIntn[dfIntn$mcBv==0 & dfIntn$targV==16,]
(p <- ggplot(iss, aes(x=trial, y=intn, colour=targEcc, linetype=stairStart)) +
facet_wrap(subj ~ targTpeak, ncol=6) + geom_point() + geom_line() + ylim(-2,0))
ss <- sumFn(thresh[thresh$maskType=='standard' & thresh$targV==0,])
ssAve <- ss[ss$subj=='average',]
ssIndiv <- ss[ss$subj!='average',]
pTargStatMaskSlow <- themefy(plotAve(ssAve)) + ylim(yLim)
(plotIndiv(ssIndiv) + ylim(-2,0))
iss <- dfIntn[dfIntn$mcBv==10 & dfIntn$targV==0,]
(p <- ggplot(iss, aes(x=trial, y=intn, colour=targEcc, linetype=stairStart)) +
facet_wrap(subj ~ targTpeak, ncol=6) + geom_point() + geom_line() + ylim(-2,0))
ss <- sumFn(thresh[thresh$maskType=='standard' & thresh$targV>0,])
ssAve <- ss[ss$subj=='average',]
ssIndiv <- ss[ss$subj!='average',]
pTargDynaMaskSlow <- themefy(plotAve(ssAve)) + ylim(yLim)
(plotIndiv(ssIndiv) + ylim(-2,0))
iss <- dfIntn[dfIntn$mcBv==10 & dfIntn$targV==16,]
(p <- ggplot(iss, aes(x=trial, y=intn, colour=targEcc, linetype=stairStart)) +
facet_wrap(subj ~ targTpeak, ncol=6) + geom_point() + geom_line() + ylim(-2,0))
ss <- sumFn(thresh[thresh$maskType=='fast' & thresh$targV==0,])
ssAve <- ss[ss$subj=='average',]
ssIndiv <- ss[ss$subj!='average',]
pTargStatMaskFast <- themefy(plotAve(ssAve)) + ylim(yLim)
(plotIndiv(ssIndiv) + ylim(-2,0))
iss <- dfIntn[dfIntn$mcBv==60 & dfIntn$targV==0,]
(p <- ggplot(iss, aes(x=trial, y=intn, colour=targEcc, linetype=stairStart)) +
facet_wrap(subj ~ targTpeak, ncol=6) + geom_point() + geom_line() + ylim(-2,0))
ss <- sumFn(thresh[thresh$maskType=='fast' & thresh$targV>0,])
ssAve <- ss[ss$subj=='average',]
ssIndiv <- ss[ss$subj!='average',]
pTargDynaMaskFast <- themefy(plotAve(ssAve)) + ylim(yLim)
(plotIndiv(ssIndiv) + ylim(-2,0))
iss <- dfIntn[dfIntn$mcBv==60 & dfIntn$targV==16,]
(p <- ggplot(iss, aes(x=trial, y=intn, colour=targEcc, linetype=stairStart)) +
facet_wrap(subj ~ targTpeak, ncol=6) + geom_point() + geom_line() + ylim(-2,0))
grid.arrange(pTargStatMaskStat + #theme(legend.position='none') +
ggtitle('a: static mask, static target'),
pTargDynaMaskStat + ggtitle('d: static mask, dynamic target'),
pTargStatMaskSlow + ggtitle('b: standard mask, static target'),
pTargDynaMaskSlow + ggtitle('e: standard mask, dynamic target'),
pTargStatMaskFast + ggtitle('c: fast mask, static target'),
pTargDynaMaskFast + ggtitle('f: fast mask, dynamic target'),
ncol=2, nrow=3) #widths=c((1-w)/2,w/2))
## Warning: Removed 3 rows containing missing values (geom_linerange).
cent <- function(v){
v <- apply(v,2,function(x){
x <- x - mean(unique(x))
x <- x / max(x)
})
}
dsc <- ds
centCols <- c('dom','targTpeak','stairStart','targEcc','maskV','targV')
dsc[,centCols] <- cent(dsc[,centCols])
# dsc$maskV <- (dsc$maskV / max(dsc$maskV)) #*2 - 1
dsc$maskV_c <- (dsc$maskV + 1) / 2
dsc$targV_c <- (dsc$targV + 1) / 2
# dsc$targEcc_c <- (dsc$targEcc + 1) / 2
dsc$targEcc_c <- round((dsc$targEcc + 1) / 2,0)
dsc$targTpeak_c <- (dsc$targTpeak + 1) / 2
# dsc$subj <- as.factor(dsc$subj)
dsc$maskType = factor(dsc$maskType, c('static','standard','fast'))
head(dsc)
## subj dom sess maskV targTpeak targEcc targV stairStart thresh cond maskType condSplit maskV_c
## 1 12 1 2 1.0000000 1 -1 1 1 -1.378333 2_20 fast 1.0000000
## 2 12 1 2 -0.6363636 1 -1 1 -1 -1.375000 2_20 static 0.1818182
## 3 12 1 2 1.0000000 -1 -1 -1 -1 -1.263333 2_20 fast 1.0000000
## 4 12 1 2 1.0000000 0 -1 1 -1 -1.341667 2_20 fast 1.0000000
## 5 12 1 2 -0.3636364 1 -1 -1 1 -1.363333 2_20 standard 0.3181818
## 6 12 1 2 -0.6363636 1 -1 -1 1 -1.636667 2_20 static 0.1818182
## targV_c targEcc_c targTpeak_c
## 1 1 0 1.0
## 2 1 0 1.0
## 3 0 0 0.0
## 4 1 0 0.5
## 5 0 0 1.0
## 6 0 0 1.0
pvalfn(lmer(thresh ~ dom + stairStart + sess + maskType * targTpeak_c * targEcc_c +
(1|subj), data=dsc[dsc$targV==-1,]))
## estm se df tVal pVal star
## (Intercept) -1.248 0.119 10.910 -10.497 0.000 ***
## dom 0.048 0.102 5.933 0.474 0.653
## stairStart 0.055 0.011 266.956 5.086 0.000 ***
## sess -0.141 0.031 271.534 -4.631 0.000 ***
## maskTypestandard 0.457 0.059 266.956 7.691 0.000 ***
## maskTypefast 0.292 0.059 266.956 4.914 0.000 ***
## targTpeak_c 0.005 0.065 266.956 0.083 0.934
## targEcc_c 0.380 0.059 266.956 6.397 0.000 ***
## maskTypestandard:targTpeak_c -0.125 0.092 266.956 -1.358 0.176
## maskTypefast:targTpeak_c -0.088 0.092 266.956 -0.954 0.341
## maskTypestandard:targEcc_c -0.213 0.084 266.956 -2.534 0.012 *
## maskTypefast:targEcc_c -0.220 0.084 266.956 -2.612 0.010 *
## targTpeak_c:targEcc_c -0.110 0.092 266.956 -1.199 0.232
## maskTypestandard:targTpeak_c:targEcc_c 0.187 0.130 266.956 1.439 0.151
## maskTypefast:targTpeak_c:targEcc_c 0.234 0.130 266.956 1.799 0.073 .
pvalfn(lmer(thresh ~ dom + stairStart + sess + maskType * targTpeak_c * targEcc_c *
targV_c + (1|subj), data=dsc))
## estm se df tVal pVal star
## (Intercept) -1.391 0.104 10.204 -13.378 0.000 ***
## dom 0.035 0.091 6.000 0.388 0.712
## stairStart 0.044 0.008 543.000 5.353 0.000 ***
## sess -0.046 0.016 543.000 -2.813 0.005 **
## maskTypestandard 0.457 0.063 543.000 7.242 0.000 ***
## maskTypefast 0.292 0.063 543.000 4.627 0.000 ***
## targTpeak_c 0.005 0.069 543.000 0.078 0.938
## targEcc_c 0.380 0.063 543.000 6.024 0.000 ***
## targV_c 0.252 0.063 543.000 3.989 0.000 ***
## maskTypestandard:targTpeak_c -0.125 0.098 543.000 -1.279 0.202
## maskTypefast:targTpeak_c -0.088 0.098 543.000 -0.899 0.369
## maskTypestandard:targEcc_c -0.213 0.089 543.000 -2.386 0.017 *
## maskTypefast:targEcc_c -0.220 0.089 543.000 -2.460 0.014 *
## targTpeak_c:targEcc_c -0.110 0.098 543.000 -1.129 0.260
## maskTypestandard:targV_c -0.140 0.089 543.000 -1.566 0.118
## maskTypefast:targV_c 0.082 0.089 543.000 0.924 0.356
## targTpeak_c:targV_c -0.017 0.098 543.000 -0.175 0.861
## targEcc_c:targV_c -0.194 0.089 543.000 -2.176 0.030 *
## maskTypestandard:targTpeak_c:targEcc_c 0.187 0.138 543.000 1.355 0.176
## maskTypefast:targTpeak_c:targEcc_c 0.234 0.138 543.000 1.694 0.091 .
## maskTypestandard:targTpeak_c:targV_c 0.090 0.138 543.000 0.650 0.516
## maskTypefast:targTpeak_c:targV_c 0.100 0.138 543.000 0.720 0.472
## maskTypestandard:targEcc_c:targV_c 0.161 0.126 543.000 1.278 0.202
## maskTypefast:targEcc_c:targV_c 0.165 0.126 543.000 1.310 0.191
## targTpeak_c:targEcc_c:targV_c 0.082 0.138 543.000 0.592 0.554
## maskTypestandard:targTpeak_c:targEcc_c:targV_c -0.136 0.196 543.000 -0.697 0.486
## maskTypefast:targTpeak_c:targEcc_c:targV_c -0.185 0.196 543.000 -0.946 0.345